import numpy as np
import matplotlib.pyplot as mp
import load_foil as lf


EPSILON = 10**-5


def splint(xs, ys, y2s, j):
  "Fonction determinant les coefficients du polynome interporlateur"
  def sp_fct(x):
    A = (xs[j+1]-x)/(xs[j+1]-xs[j])
    B = (x-xs[j])/(xs[j+1]-xs[j])
    C = (((A**3) - A) * ((xs[j+1] - xs[j])**2))/6
    D = (((B**3) - B) * ((xs[j+1] - xs[j])**2))/6

    y = A*ys[j] + B*ys[j+1] + C*y2s[j] + D*y2s[j+1]
    return y

  d = ((xs[j+1]**2)*xs[j]*(2*y2s[j]+y2s[j+1]) + xs[j+1]*(6*ys[j]-(xs[j]**2)*(y2s[j]+2*y2s[j+1])) - 6*xs[j]*ys[j+1]) / (6*(xs[j+1]-xs[j]))
  c = (-2*(xs[j+1]**2)*y2s[j] - (xs[j+1]**2)*y2s[j+1] - 2*xs[j+1]*xs[j]*y2s[j] + 2*xs[j+1]*xs[j]*y2s[j+1] + (xs[j]**2)*y2s[j] + 2*(xs[j]**2)*y2s[j+1] - 6*ys[j] + 6*ys[j+1]) / (6*(xs[j+1]-xs[j]))
  b = (xs[j+1]*y2s[j] - xs[j]*y2s[j+1]) / (2*xs[j+1] - 2*xs[j])
  a = (y2s[j+1]-y2s[j]) / (6*(xs[j+1]-xs[j]))
  poly_coef = np.array([a,b,c,d])

  return np.array([sp_fct,poly_coef])

def spline_lambda(xs, ys, y2s):
  "Fonction retournant un tableaux contenant des lambda et les coefficients associees pour chacun des polynomes sur l'ensemble de l'intervale de definition"
  n = np.shape(xs)[0] 

  spline_lambda = []
  for i in range(0,n-1):
    spline_lambda.append(splint(xs, ys, y2s, i))
  return np.array(spline_lambda)  


def second_derivative(x, y):
  "Fonction calculant les derivees secondes"
  n = np.shape(x)[0]
  u = np.zeros(n)
  y2 = np.zeros(n)
  
  y2[0] = 0.0
  u[0] = 0.0
  

  for i in range(1, n-1):
    sig = (x[i]-x[i-1])/(x[i+1]-x[i-1])
    p = sig * y2[i-1] + 2.0
    y2[i] = (sig - 1.0)/p
    u[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i] - y[i-1])/(x[i]-x[i-1])
    u[i] = (6.0*u[i])/(x[i+1]-x[i-1]) - (sig*u[i-1])/p

  
  qn = 0.0
  un = 0.0

  y2[n-1] = (un-qn*u[n-2])/(qn*y2[n-2]+1.0)
  for k in np.arange(n-2, 0, -1):
    y2[k] = y2[k]*y2[k+1] + u[k]
  return y2
  

def spline(x, y):
  "Fonction interpolant un tableau de point x,y par spline cubique"
  y2s = second_derivative(x, y)
  return spline_lambda(x, y, y2s)


def plot_spline():
  "Fonction tracant la fonction interpolee"
  (ex,ey,ix,iy) = lf.load_foil("naca747a315_short.dat")

  e_spline_value = spline(ex, ey)
  for i in range(0, len(ex) - 1):
    mp.plot(np.arange(ex[i], ex[i+1], 0.0001), e_spline_value[i][0](np.arange(ex[i], ex[i+1], 0.0001)))

  i_spline_value = spline(ix, iy)
  for i in range(0, len(ix) - 1):
    mp.plot(np.arange(ix[i], ix[i+1], 0.0001), i_spline_value[i][0](np.arange(ix[i], ix[i+1], 0.0001)))


  #mp.plot(ex, ey)
  #mp.plot(ix, iy)
  mp.ylabel('y')
  mp.xlabel('x')
  mp.ylim(bottom=-0.2, top= 0.2)
  mp.show()
  mp.clf()


if __name__ == "__main__":
  plot_spline()
